Exact Mapping Noisy van der Pol Type Oscillator onto Quasi-symplectic Dynamics 
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We find exact mappings for a class of limit cycle systems with noise onto quasi-symplectic dy- 
namics, including a van der Pol type oscillator. A dual role potential function is obtained as a 
component of the quasi-symplectic dynamics. Based on a stochastic interpretation different from 
the traditional Ito's and Stratonovich's, we show the corresponding steady state distribution is the 
familiar Boltzmann-Gibbs type for arbitrary noise strength. The result provides a new angle for 
understanding processes without detailed balance and can be verified by experiments. 
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Noise disturbed limit cycle dynamics is now attracting 
considerable attention in the physics community [1-3]. A 
direct reason is that ubiquitous real systems can be mod- 
eled by them, e.g., from cell cycle [4], circadian phenom- 
ena [5] to chemical reaction [6] and oscillating electrical 
circuit [7] . The dynamics itself is a touchstone to study 
nonlinear dissipative process in the absence of detailed 
balance. Due to the difficulty arising out of nonlincarity 
and stochasticity, approximated methods based on phase 
reduction and weak noise perturbation are proposed from 
former studies [1, 2], but exact results are rarely seen in 
literature. In this Letter, we examine noised limit cy- 
cles through exactly mapping them onto quasi-symplectic 
[8] dynamics (defined in Eq. (3)) valid for an arbitrary 
noise strength. A potential function is obtained therein 
as a component of the quasi-symplectic dynamics. These 
functions serve as both a Lyapunov function of the de- 
terministic part dynamics [9] and a Hamiltonian leading 
to Boltzmann-Gibbs type steady state distribution of the 
stochastic process [10]. 

The existence of such potential function for processes 
without detailed balance is still suspected [3, 11], a spe- 
cific argument is that: for a limit cycle system with non- 
constant velocity along the cycle, the dual role potential 
does not exist. We begin with such an example. The van 
der Pol oscillator [7] is a famous limit cycle dynamics, 
here we examine a stochastic version with a multiplica- 
tive noise ({c[,t) = (Ci(q, t), (2(1, t)) T (the superscript r 
denotes the transpose of a matrix) and a higher order 
term h(q\); 
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When h(qi) = 0, the deterministic part of the dynamics 
reduces to the van der Pol oscillator. A specific system 
we would like to illustrate is h{q\) = ^ 2 q\/A — [i 2 q\/\6 
[12], we can observe from Fig. 1, the deterministic dy- 
namical behavior of the system. It has a rotationally 
non-symmetric limit cycle and a position dependent ve- 
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FIG. 1. Upper panel: Trajectories (deterministic dynamics) 
for the system Eq. (1) with ft(gi) = ^qf/A-^ql/lG (jj, = 1). 
Lower panel: Comparison of two systems, the van der Pol 
Oscillator h(qi) — is in blue, the system in the upper panel 
is in red. 



locity along the cycle. We propose three questions here: 

(1) . Whether a potential function (Lyapunov function) 

can be explicitly constructed for the system in the 
upper panel of Fig. 1? 

(2) . Under which kind of stochastic integration can the 

SDEs in Eq. (1) have a Boltzmann-Gibbs steady 
state distribution for an arbitrary noise strength? 

(3) . What is the significance of constructing such poten- 

tial functions? 
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To answer these questions, we should first briefly re- 
view a framework [10]. Consider the Langevin equation 
[13]: 



q = f (q) + N{q)Z(t) , 



(2) 



where q, f are n-dimensional vectors and f is a non- 
linear function of the state variable q. The noise term 
is fc-dimcnsional Gaussian white with the zero mean, 
= 0, and the covariance (£(t)£ T (i')) = S(t - t')I k . 
The notation S(t — t') is the Dirac delta function, (...) 
indicates the average over noise distribution, is the k- 
dimcnsional identity matrix. The element of the n x k 
matrix iV(q) can be a nonlinear function of q, then the 
noise considered in this framework can be a general mul- 
tiplicative noise. This matrix is further described by: 
N(q)N T '(q) = 2eD(q), the constant e quantifying the 
noise strength and D(q) is a n x n positive semi-definite 
diffusion matrix. Note that the noise may have less than 
n independent components k < n, leading to the zero 
eigenvalue of D(q). During the study of a biological 
switch [10], a quasi-symplectic dynamics equivalent to 
Eq. (2) was discovered: 



[5(q) + A(q)]q = -V0(q) + #(q)£(i) 



(3) 



The term S(q) is a positive semi-definc matrix, — 5(q)q 
denotes a frictional force; the term A(q) is an anti- 
symmetric matrix representing an embedded symplectic 
structure, — A(q)q is a rewritten of the Lorentz force 
eq x B in 2 or 3 dimensional cases, and also a gener- 
alization to higher dimensions; the scalar function </>(q) 
is a potential function, e.g., the electrostatic potential, 
lying at the core of the discussion in this Letter. The 
matrix N(q) is constrained by the fluctuation-dissipation 
theorem [14]: N(q)N T (q) = 2eS(q). 

A corresponding Fokker-Planck equation for (3) (there- 
fore for Eq. (2)) can be obtained with physical signifi- 
cance (a zero mass limit) [15]: 

d t p(q, t) = V ■ [D(q) + Q(q)] • [V0(q) + eV]p(q, t) , (4) 

where V in V0(q) does not operate on p(q, t); D(q) 
is the diffusion matrix; the matrix Q(q) is anti- 
symmetric and can be calculated from the relation 
[S{q) + A(q)] [D(q) + Q(q)] = I n . Equation (4) has the 
Boltzmann-Gibbs distribution with the potential <f>{q) as 
a steady state solution: 



p(q,t -> oo) 



1 



■ exp 



(q) 



(5) 



where Z s = J rfqexp {cf>(q)/e} is the partition func- 
tion. The 1-d case has been verified by an experi- 
ment [16]. The probability current density j(q, t) — 
(ji(q, t), . . . , j„(q, t)) T is commonly defined as: 

ji(q, t) = /i(q)p(q, t) - dj [eDy (q)/>(q, *)] (6) 



where /j(q) = fi(q)+ed j [D i j(q)+Q il (q)}, fi(q) is the zth 
component of the vector valued function f(q) in Eq. (2), 
Dij(q) and Qij(q) are the elements of the matrices D(q) 
and Q(q) in Eq. (4). In steady state, the probability dis- 
tribution is given by Eq. (5). We have V ■ j s (q) = 0, but 
j s (q) is usually not zero. One can check that Q — is a 
sufficient condition for j s = 0; but when Q(q) ^ then 
generally j s (q) 7^ 0, since dj [Q»j(q)p a (q)] 7^ 0. There- 
fore the framework encompasses the cases without de- 
tailed balance. The term "detailed balance" means the 
net current between any two states in the phase space is 
zero [17], identical to that for Markov process in mathe- 
matics. The dynamics studied in this Letter corresponds 
to the non-detailed balance cases discussed in [18] as well. 

Equation (4) then defines a stochastic interpretation 
for the Langevin equation Eq. (2), that is the A- type 
integration. The relation between friction S(q) and dif- 
fusion -D(q) in the absence of detailed balance condition 
is named as the generalized Einstein relation in [10], in 
1-d case, it reduces to the well-known Einstein relation. 
A detailed derivation of Eq. (4) has been provided in [15] . 

The Langevin equation (2) can also be considered as 
a composition of a deterministic dynamics q = f (q) and 
a fluctuation Af(q)£(t). For the deterministic dynamics, 
the time derivative of the potential function <f>(q) along 
a trajectory is: 

#(q) 



dt 



= • f (q) = -V0(q) • [D(q) + Q(q)] • V0(q) 
= -V<Xq) • D(q) ■ V<Kq) < , 



since D(q) is nonncgative and symmetric. It means that 
the potential along the trajectory is non-increasing and 
has the local extreme values at fixed points, limit cycles 
or more complex attractors. Hence, potential function 
<ft{q) servers as a Lyapunov function [9] for the determin- 
istic dynamics q = f(q). When noise exists, the system 
can deviate from the attractors, transferring from one 
stable state to another, once the system enters a new 
basin of attraction, the deterministic force (dissipation) 
will drive the system to the corresponding attractor. It 
shows that noise is sometimes essential in mathemati- 
cal modeling of reality and plays a driving role in phase 
transitions. 

For all planar rotationally symmetric limit cycle dy- 
n amics (r epresented below in polar coordinate (q = 
\J q\ + 92 1 0) f° r simplicity) with the diffusion matrix 
D = D l2 (Dq is the constant diffusion coefficient): 

q = R(q)+Cg(q,e,t) 
9 = ij(q)+Ce(q,0,t) ' 



(7) 



we can provide the exact construction for the potential 
function (note that the following result should be trans- 
formed back to Cartesian coordinate): 



^(q) = - 



R(q)dq 



(8) 
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and corresponding dynamical components: 



x x y = (xiUj — Xjyi) nxn , the result is an n x n matrix): 



S(q) 
Q(q) 



R{qf 



D [R(q) 2 + q 2 ^(q) 2 
D [R(q) 2 + q 2 ^(q) 2 

g^(g)A) ( o l 

fl(g) ^ -1 



1 
1 

1 

-1 



(9) 



Note that for the weak noise limit, when Dq — > 0, this 
construction is still valid for the deterministic dynamics 
by an arbitrarily setting D (see also [19]). 

More generally, we can extend our method to rota- 
tionally non-symmetric systems with multiplicative noise 
through coordinate transformations (reversible continu- 
ous mapping, can be nonlinear) a of the deterministic 
dynamics: 



qi = /i(q) 
92 = .Mq) 



< u = h(u,v) ^ 
v = h{u,v) 



q = p{q)P{q,6) 
= <p(q,6) 



where P(q 1 9) 0. The potential function is given by 

<K?l,«2) < (U ' V) 4>{u,v) ( P ° lar (q) (p p (q) = - p(q)dq 



since d<t>{qi , q 2 ) / 'dt = [d(j> p {q) / dq] q = -p 2 (q)P(q,9) < 0. 
The protocol of the construction: First, for the de- 
terministic dynamics (91,(72), find out the transforma- 
tion a : (91,92) — > {u, v), calculate the dynamics under 
(u,v), that is the u — fi(u,v) and v — f2(u,v); Sec- 
ond, rewrite the obtained dynamics in polar coordinates 
(u, v) — > (q,9), if the dynamics can be expressed as the 
requested form above, the potential function can be con- 
structed as 4> p (q); Third, transform (j> p {q) back to 4>(u, v), 
and finally to (j>{q\, 92). 

Once potential function (Lyapunov function) </>(q) for 
the deterministic dynamics q = f(q) is obtained, there 
are different ways to construct the dynamical compo- 
nents, one particular setting is provided in [19] (the bi- 
nary operator of two n dimensional vectors is defined as 



5(q) = -^/, A(q) 



Q(q) 



f • f 
ff 



V0 x f 

ff 

2 



(V0 x f) 



V0-f (V<£ • f ) (V<f> • V<£) 
V</> x f 



(10) 



We should note that although the potential function 4> 
is invariant during coordinate transformation, the dy- 
namical components, the matrices S, A, D, Q vary in 
different coordinates, but a straightforward formulation 
can be achieved by multiplying Jacobian matrix of the 
transformation. 

Back to the example in Eq. (1) with (Ci(q, t), C2(q, t)Y 
= iV(q) • (£i(t),&(t)) T and h( qi ) = p 2 qf/A - p 2 q\/l&. 
The deterministic part is a Lienard equation similar to 
the famous van der Pol Oscillator (0 < p < 2, sec 
Fig. 2) [12]. Through a nonlinear coordinate transfor- 
mation <j~ l : u = 91 and v = 92 — pqi + M9i/4, we obtain 
the dynamical system: 



f u = pu - f u 3 + v ^ f 



S(4-9 



2 ^ 9 cos 2 



-1 — p cos $ sin 9 



with ,0(9) = (4 - q 2 )q and P(g, 0) = p cos 2 6>/4 > (p > 
0). Therefore, we can provide an exact construction of 
potential function for Eq. (1) (sec Fig. 2): 



0(q) 



a* 3 

T9i 



9? + (92 - imi + 

?i + (<fc-Wi + ^?i) -8 



(11) 



We note that the potential function Eq. (11) has the min- 
imal value at the stable limit cycle 91 = pq\ — jqf ± 
y/ A — q\ and a local maximum value at the unstable fixed 
point (0, 0). Expressions for other dynamical components 
can be provided through Eq. (10). We use the represen- 
tation with (u = 91, v = 92 — j«9i + M9i/4) and J(q) the 
Jacobian matrix d(u,v)/d(qi,q2): 



m = " (4 4 ,g;;:; v -""> Tj w - - - (, "'it';' tH * ( ° i ) ^ • 

I 



The result obtained can be understand as: The stochas- 
tic dynamics Eq. (1) with the position dependent dif- 
fusion matrix Z?(q) given in Eq. (12) has the explic- 



itly constructed potential function </>(q) (Eq. (11)) and 
a corresponding Boltzmann-Gibbs steady state distribu- 
tion (Eq. (5)). For the matrix <2(q), one can check 
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FIG. 2. Potential function Eq. (11) with fi = 1: The color 
blue denotes a lower potential value, and the red means higher 
value. The graph is drawn below a preset upper bound value 
1, the phase variables are qi and q-z. 

9j [Qij(<l)p s (q)] =^ 0, leading to the absence of de- 
tailed balance. The stochastic integration used is the 
A-type (see Eq. (4)) different from traditional Ito's or 
Stratonovich's [15]. Note that when approaching to the 
limit cycle (4 — u 2 — v 2 ) — > 0, the force induced by the 
potential gradient goes to zero; the Lorentz force matrix 
A(q) goes zero in the same order and changes its sign 
at the limit cycle (since < /i < 2); the friction matrix 
S(q) goes to zero in a higher order. The dynamics at the 
limit cycle is no longer dissipative but conserved in this 
limit, reaching a stable cycle. Thus the potential value 
should be equal on limit cycles where the system is con- 
served. We note that this is consistent with the definition 
of a Lyapunov function [9], since the particle is moving 
repeatedly along the cycle, the same as the conserved 
system moving along the Hamiltonian. The singularity 
problem for this construction has been discussed in [19]. 
Previous works focus more on the diffusion matrix, ig- 
noring the important role played by the friction matrix 
S(q) and the Lorentz force matrix A(q). 

The significance of constructing potential functions de- 
fined in this framework can be viewed from two aspects: 
First, there is a simple and direct correspondence be- 
tween stochastic and determinist dynamics for arbitrary 
noise strength (not only under weak noise limit). For ex- 
ample, attractors like fixed points or limit cycles do not 
change their position in the phase space. This enables a 
straightforward use of the dynamical analysis for the de- 
terministic dynamics in the presence of noise. Therefore, 
the calculation of the transition probability from one sta- 
ble fixed point qi to another one through a saddle point 
q2 is generally formulated as cx exp [ — |<^(qi) — <^(q 2 )|/e]; 
Second, potential function obtained here serves also a 



Lyapunov function for the deterministic dynamics. The 
framework then contributes possible new approaches for 
the largely unsolved problem in engineering: construct- 
ing Lyapunov function for general nonlinear dynamics. 
In addition, the A-type stochastic integration can be ap- 
plied directly in the study of phase reduction. For con- 
ventional phase reduction method [1] , A-type integration 
does not lead to the noise-induced frequency shift (NIFS). 

In conclusion, we have exactly mapped a class of limit 
cycle systems with noise onto quasi-symplectic dynamics, 
from rotationally symmetric to non-symmetric systems, 
a specific example is a van der Pol type oscillator. These 
systems have been analyzed through the dynamical com- 
ponents S(q), A(q) and </>(q) obtained by the mapping. 
Using A-type integration, the steady state distribution of 
these systems is Boltzmann-Gibbs type, consistent with 
the knowledge from statistical physics. Since new mea- 
suring techniques for Brownian motion is available (e.g. 
[16]), the theoretical results here may be verified by ex- 
periments directly in the near future. 
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